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Modelling the Bioclimatic Niche of a Cohort of Selected Mite Species (Acari, Acariformes) Associated 
with the Infestation of Stored Products. Tytar, V. M., Oksentyuk, Ya. R. — In this study an attempt 
is made to highlight important variables shaping the current bioclimatic niche of a number of mite 
species associated with the infestation of stored products by employing a species distribution modeling 
(SDM) approach. Using the ENVIREM dataset of bioclimatic variables, performance of the most robust 
models was mostly influenced by: 1) indices based on potential evapotranspiration, which characterize 
ambient energy and are mostly correlated with temperature variables, moisture regimes, and 2) strong 
fluctuations in temperature reflecting the severity of climate and/or extreme weather events. Although 
the considered mite species occupy man-made ecosystems, they remain more or less affected by the 
surrounding bioclimatic environment and therefore could be subjected to contemporary climate change. 
In this respect investigations are needed to see how this will affect future management targets concerning 
the safety of food storages. 
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Introduction 


More than 50 species of mites have been found in association with stored products (Clemmons, Taylor, 
2016). Mites are important pests in stored grain commodities not only due to the damage caused by their feed- 
ing, but their presence may lead to rejections resulting in a loss in revenue. In addition they are highly allergenic 
and have the potential to vector mycotoxin producing fungi within a storage environment which may result in 
associated health risks if stores are not properly managed (Collins, 2012). Stored-grain mites are, for the most 
part, cosmopolitan (Freeman, 1973). Most stored-grain arthropods are thought to have originated in India or 
the Middle East and have been transported throughout the world by human movement of foods (Cotton, 1956). 
The main limitations to the spread of stored-grain pests, besides the export and import of goods, and the regula- 
tory efforts of importing countries, are temperature (T) and relative humidity (RH) (i. e. climate) (White et al., 
2011). Mites are very common pests of stored grain and oilseeds, particularly in temperate areas with high RH, 
although they can also infest stored products in relatively dry and warm climates (Palyvos et al., 2008). 

The most important factors that influence mite population development appear to be combinations of T, 
RH and moisture content (MC). A rise in T, RH and/or MC increases mite productivity, whereas a decrease 
in T, RH and/or MC content decreases productivity, with interactions occurring between the variables (Dunn, 
2003). 
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Different mite species also have different requirements for development, although there is generally no 
upper limit for humidity (Cunnington, 1969). Acarus siro Linnaeus, 1758 has a greater rate of development in 
grain of higher MC and is less frequent in drier warmer climates (Cunnington, 1976). Under laboratory condi- 
tions, the lower and upper T limits for A. siro are around 2.5 and 31 °C respectively; with a lower RH limit of 
62.5 % (dependent on T) (Cunnington, 1965), although Aspaly et al. (2007) found an upper T limit of 38 °C 
(at 85 % RH). Lepidoglyphus destructor (Schrank, 1781) is less tolerant of low humidities than A. siro and can- 
not complete development below much 65 % RH (dependent on T) (Cunnington, 1976). Its lower and upper 
T limits, are around 3 and 34 °C respectively (Cunnington 1976). Tyrophagus putrescentiae (Schrank, 1781) is 
more tolerant of high temperatures and can survive and breed readily above 30 °C; but is less tolerant of low 
temperatures and cannot develop much below 10 °C (Cunnington, 1969). As with the other storage mites, 
T. putrescentiae develops most rapidly at conditions of high RH, with a lower limit of 65 % (dependent on T) 
(Cunnington, 1969). The optimum conditions for rapid development in the laboratory are 25 °C and 90 % RH 
for A. siro and L. destructor, and 30 °C and 90 % RH for T. putrescentiae (Cunnington, 1976). Due to the wide 
range of preferred temperatures, A. siro and L. destructor can be classified as eurythermic, whereas T. putres- 
centiae is stenothermic and more thermophilous than the other species (Hubert et al., 2010). The species T. pu- 
trescentiae obviously gravitates in its life strategy to an r-strategy. It is thermophilic, is first to appear in stored 
products, rapidly increases in numbers, quickly spreads over the entire substrate, tends to consume substrates 
rich in nutrients, because of peculiarities of its digestive enzymes, and is first to disappear (Akimov, 1985). In 
addition, this species has no hypopial stage and is oriented towards establishing itself in new substrates. On the 
other hand, A.M. Cunnington considers limits of T and RH within which A. siro is able to complete its deve- 
lopment are fairly restricted and this explains why, although often loosely described as a cosmopolitan species, 
its occurrence and distribution is largely confined to countries with a cool, moist climate (Cunnington, 1965). 

Mite species can form a hypopial stage, which does not feed and can be resistant to unfavourable condi- 
tions such as desiccation (Griffiths, 1964) and low temperatures. For instance, there is a report of their survival 
at temperatures of -18 °C (Sinha, 1964). Interestingly, humidity is the main factor in bringing about termina- 
tion of the hypopial stage which increases as levels rise above 70 % (Stratil, Knulle, 1985). 

Mite growth is affected by a combination of physical and biological factors which have been established 
under current climatic conditions. These factors interact in a complex manner (Sinha, Wallace, 1973) and 
knowledge of the factors important in regulating mite development and infestations can allow for the most ap- 
propriate control measures to be undertaken (Collins, 2012). However, global climate patterns have shown to 
change notably (IPCC, 2018) and grain storage itself can be heavily affected by climate change. The conclusion 
that seasonal variations would have an effect on mite population (Sinha, 1968) indicates that climate change 
would also have similar effects (Moses et al., 2015). 

Because of these concerns there is a need in developing approaches aiming to understand the relative im- 
portance of current factors as a baseline in shaping mite communities. One tool to do this is species distribution 
modeling (SDM) which employs suitability indices. Suitability indices describe the relationship between habitat 
suitability score and a given environmental variable of a target species. Habitat suitability is a way to predict the 
suitability of habitat at a certain location for a given species or group of species based on their observed affinity 
for particular environmental conditions (Yi et al., 2016; Ma, Sun, 2018). 

SDMs are empirical tools in ecology, biogeography, natural resource management, and ecosystem man- 
agement (Franklin, 2009). Among various SDMs, Maxent is highly regarded because it computes the probabi- 
lity distribution of a species by using maximum entropy rules (the greatest uniform distribution) and statistical 
mechanics (Elith et al., 2006). The model has the advantage that it needs only environmental information and 
species occurrences data (Elith et al., 2011). MaxEnt is most important for modelling the suitability distribution 
of limited occurrence data and predictions may be modelled using a very few sample occurrences (Hernandez 
et al., 2008); sample sizes as low as n = 10 with Maxent performed well (Wisz et al., 2008; Kumar, Stohlgren, 
2009). 

In this study we make an attempt through building reasonable and correct (according to evaluation re- 
sults) SDMs to highlight the important variables shaping the current bioclimatic niche of selected mite species. 

It is common that SDMs employ climatic variables, however climatic variables derived from climate da- 
tasets at coarse scale may not appropriately account for anthropogenic influences on microclimate (Feder- 
man et al., 2013; Varner, Dearing, 2014), particularly in human-made structures, which may be buffered in 
terms of temperature and moisture from the outdoor environment. For instance, this may apply to an ar- 
ray of mites associated with grain storage facilities, barns, houses, even apiaries. Nevertheless, grain mites are 
more common in temperate regions with cool moist climates (Palyvos et al., 2008), showing the significance of 
the macroclimate. In terms of the macroclimate, climates which experience mild winters and high RH offer a 
particular challenge in protecting commodities from mite infestation. Even within the storage facility surface 
temperatures generally mimic the ambient conditions with moisture changes at the surface lagging about 4h 
behind the ambient RH (Sinha, 1973; Burrell, 1979 as cited in Armitage, Cook, 1999). There are considerations 
that differences in the composition of house dust mite species and their abundances may be attributed to the 
geo-climatic conditions that prevail in a given area (Lang, Mulla, 1977; Mumcuoglu et al., 1999). Studies have 
shown that honeybee pests (for example, Varroa destructor Anderson et Trueman, 2000) can survive only in 
certain optimal bioclimatic conditions. For instance, the optimal temperature, humidity, precipitation, altitude 
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and biomass/net primary productivity ranges for different honeybee pests can vary significantly (Makori et al., 
2017); the quoted authors concluded that honeybee pests could be modelled using coarse scale bioclimatic data 
which were most relevant in all model results. Nonetheless, we have to accept the notion that microclimates 
are supposed to modulate the responses to the macroclimate, but currently the extent of such modulation in a 
coarse scale modelling exercise is difficult to verify. 


Material and methods 


Study area 

The area is geographically located between 49°26’ to 52°1' N latitudes and 25°57’ to 29°43’ E longitudes, 
with an altitude ranging up to 319 m. 

Using the Koppen-Geiger climate classification (McMahon, 2007), the type of climate in the study area is 
classified as “humid continental” and represents its “cool summer” version. It has little warming or precipitation 
effects from the northern Atlantic. The cool summer subtype is marked by mild summers, long cold winters and 
less precipitation than the hot summer subtype, however, short periods of extreme heat are not uncommon. 

The historic average annual temperature in the study area is 7.17 °C and the average annual precipitation — 
636 mm. Average monthly temperature and precipitation, as well as average monthly relative humidity, 
substantially vary across the year (fig. 1-3). July on average is the warmest month (18.7 °C), January — the 
coldest (-5.8 °C); July on average is the wettest month (92.2 mm), February and March are the driest (31.3 mm); 
on average in November relative humidity is the highest (75.7 %), in May it dips to 51.6 %. The average monthly 
temperature and precipitation are highly correlated (r = 0.86, p-value < 0.05), whereas the average monthly 
relative humidity is inversely correlated with both the average monthly temperature and precipitation (r = 
-0.79 and -0.61, respectively, p-value < 0.05). 

In terms of the Holdridge life zones system used to understand biome characteristics (Holdridge, 1947), 
the target area is recognized as “cool temperate wet forest”, containing pure or mixed stands of broad-leaved 
deciduous or needle-leaved evergreen tree (predominantly Scots pine, Pinus silvestris Linnaeus) growth forms, 
with a seasonal green understory of herbs. 

According to the GLC2000 landcover classification (Bartholomé, Belward, 2005), most of the study area 
is considered to consist of “Cultivated and managed areas” (36 %) and in summary around 35 % of the area is 
under various types of tree cover and/or other natural vegetation. 

The anthropogenic impact on the study area can be assessed by the Global Human Footprint Index 
(Wildlife Conservation Society, 2005), where a value of zero represents the least impacted part of the biome 
(“wildest”) and a value of 100 represents the most influenced (least “wild”) part of the biome. In our case this 
index consists of an average score of around 34, meaning the overall human pressure on the environment is 
fairly low. 
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Fig.1. Average monthly temperature. 
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Fig. 3. Average monthly relative humidity. 


Species occurrence data and climatic predictors 

The study used presence-only species collection techniques for data collection. The presence-only data 
collection technique is favoured over the presence-absence data collection technique as the latter leads to 
pseudo-absence data of species occurrence (Gormley et al., 2011). The species occurrence data used in this 
study were obtained from the various field surveys 200 a total of 63 occurrence points of 30 acaridid species 
were identified and included in this study (table 1). For modelling purposes species with occurrences = 10 
were considered. 

The species distribution points were plotted in a GIS environment using the SAGA GIS software (Conrad 
et al., 2015). 

There are many variables (including hydrological-thermal) that are thought to be relevant to species’ 
ecology and geographic distribution. In this study we used two sets of such variables. In the first place we used 
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Table 1. Occurrence of acaridid mites 









































Occur- Feeding substrate 
Family Species ee i) alee) le 
Suidasidae Suidasia nesbitti Hughes, 1948 1 + 
Acarus siro Linnaeus, 1758 35 t+ + + + + + 
A. farris (Oudemans, 1905) 10 + + + 
A. tyrophagoides (Zachvatkin, 1941) 5 + + + 
Mycetoglyphus fungivorus Oudemans, 1932 1 + + 
Tyrolichus casei Oudemans, 1910 5 + + 
Tyrophagus putrescentiae (Schrank, 1781) 16 + + to + 
T. molitor Zachvatkin, 1941 12 + 4 + + +4 
T. perniciosus Zachvatkin, 1941 8 + + + + 
T. humerosus (Oudemans, 1923) 3 + + 
Tandé T. longior (Gervais, 1844) 2 + + 
T. formicetorum Volgin, 1948 1 + 
T. mixtus Volgin, 1948 1 + 
Schwiebea nova (Oudemans, 1906) 1 + 
Neoacotyledon sokolovi (Zachvatkin, 1940) 7 + + 
Sancassania berlesei (Michael, 1903) 6 + 
S. sphaerogaster (Zachvatkin, 1937) 7 + 
S. rodionovi (Zachvatkin, 1935) 3 + 
S. mycophagus (Megnin, 1874) 1 + 
S. oudemansi (Zachvatkin, 1937) 1 + 
Rhizoglyphus echinopus (Fumouze et Robin, 1868) 3 + 
Glycyphagus domesticus (De Geer, 1778) 57 + + + + + + 
Lepidoglyphus destructor (Schrank, 1781) 55 + + + + 4 + 
L. fustifer (Oudemans, 1903) 8 + + + 
Glycyphagidae L. burchanensis Oudemans, 1903 3 + + 
L. michaeli Oudemans, 1903 3 + + 
L. pilosus (Oudemans, 1906) 2 + 
Gohieria fusca (Oudemans, 1902) 3 + + 
Chortoglyphidae Chortoglyphus arcuatus (Troupeau, 1879) 1 + 
Aeroglyphidae Aeroglyphus peregrinans (Berlese, 1892) 5 + 





Note. 1 — oilseed crops; 2 — grain; 3 — fodder; 4 — hay and straw; 5 — litter, dead bees and ambrosia 
from beehive bottoms; 6 — spoiled vegetable crops (root and tuberous crops). 


the widely accepted bioclimatic potential predictor variables for species distribution and suitability analysis 
(Hijmans et al., 2005). These bioclimatic predictors are ecologically more sensitive to differentiate the physio- 
ecological tolerances of a habitat (Thompson et al., 2009) than simple temperature and precipitation predictors 
(Graham, Hijmans, 2006; Kumar et al., 2009). Information on the bioclimatic parameters was collected as raster 
layers from the WorldClim website (http://www.worldclim.org/current) with a spatial resolution of 30 arc sec- 
onds (this is about 1 km at the equator). The dataset was prepared by the interpolation of the historical records 
(1950-2000) of monthly precipitation and monthly temperature (Hijmans et al., 2005). These 19 bioclimatic 
variables indicate a general trend of precipitation and temperature, extremity and seasonality of temperature. 
Bioclimatic predictors were used to characterize the species’ bioclimatic profile. 

Secondly, we have used a recently reconsidered (in terms of biological significance) set of 16 climatic and 
two topographic variables (the ENVIREM dataset, downloaded from http://envirem.github.io), many of which 
are likely to have direct relevance to ecological or physiological processes determining species distributions 
(Title, Bemmels, 2018). Considering the relative evenness of the terrain in the study area, topographic variables 
were excluded from the analysis. 

The most important sets of variables assessed in performance of models were identified using the Maxent 
Variable Selection package in R (developed by Jueterbock et al., 2016). The program first tests for Pearson correlation 
(set for r < 0.8), then through an iterative process tests different regularization values, correlation and contribution of 
variables (set for > 5 %) to determine the best possible set of variables by sample size corrected Akaike information 
criteria (Akaike, 1975) and AUC (see 2.3). In Maxent, the regularization parameter (B) prevents overfitting, that is, 
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Fig. 4. Response of A. siro to PETcoldQ. 


the model loses the ability of generalizing beyond the fitting data. The best possible sets of variables were selected for 
models which produced the highest AUC for the test dataset (Warren, Seifert, 2011). 


Predictive niche modelling 

The study used Maxent (Maximum Entropy, version 3.3.3k downloaded from http://www.cs.princeton. 
edu). Pseudo-absence points (used instead of true absences) were randomly generated within a convex hull 
encompassing 63 presence points. Within the settings window, a bootstrap replicate run (i. e., sampling with 
replacement) type was selected for 50 replicates with a random test percentage of 30 percent used. The option 
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Fig. 5. Response of L. destructor to PETcoldQ. 
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Fig. 6. Response of GI. domesticus to PETcoldQ. 


to use a random seed was chosen so that each replicate run would start with a random seed to ensure that a 
separate test/train dataset was used for each of the replicate models. An independent dataset was not available 
for model assessment and we did not want to partition the dataset into test and training data, and lose valuable 
training data. As an alternative, bootstrap replication splits the dataset multiple times, and in each case, 
predictive performance is assessed against the test dataset. This allows for all occurrence data to be used with a 
random partition performed with each Maxent replicate run. 

Maxent provides output data in various formats. The logistic format (values range from 0 to 1) is 
recommended, because it allows for an easier and potentially more accurate interpretation of habitat suitability 
compared to the other approaches (Baldwin, 2009). 
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Fig. 7. Response of A. siro to aridityIndexThornthwaite. 
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Fig. 8. Response of L. destructor to aridityIndex Thornthwaite. 


Maxent modeling can determine the importance of environmental variables. In one option the 
contribution for each variable can be determined by randomly permuting the values of that variable and 
measuring the resulting model performance. Better performance means that the model depends heavily on that 
variable (for explicitness values are normalized to give percentages). The second option uses a jackknife test 
and the regularized gain change during each iteration of the training algorithm. The environmental variable 
with the highest gain is considered to have the most useful information by itself, whereas the variable causing 
the largest decrease in the model’s gain contains the most information not found in the other environmental 
variables. These options are used to determine the importance of the environmental variables and to observe 
whether or not the predictive distribution responded to each variable as expected. 
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Fig. 9. Response of GI. domesticus to aridityIndexThornthwaite. 


Modelling the Bioclimatic Niche of a Cohort of Selected Mite Species... 407 





= = 
~ 


> 
in 


Probability of presence 


0.4 








+r t T T t 
23.2 23.4 23.6 23.8 24 24.2 24.4 24.6 24.8 
continentality 


Fig. 10. Response of L. destructor to continentality. 


Maxent also allows the construction of response curves to illustrate the effect of selected variables on 
habitat suitability (consequently, on the probability of occurrence and giving an idea of where for each variable, 
under the constraints and conditions of the modelling situation, the focal species has it’s optimum). These 
response curves consist of the specific environmental variable as the x-axis and, on the y-axis, the predicted 
probability of suitable conditions as defined by the logistic output. Upward trends for variables indicate a 
positive relationship; downward movements represent a negative relationship; and the magnitude of these 
movements indicates the strength of the relationship (Baldwin, 2009). 

An important part of determining the ability of niche models to predict the distribution of a species is 
having a measure of fit. The performance of the Maxent model is usually evaluated by the threshold-independent 
receiver operating characteristic (ROC) approach (calculating the area under the ROC curve (AUC) as a 
measure of prediction success). The ROC curve is a graphical method that represents the relationship between 
the false-positive fraction (one minus the specificity) and the sensitivity for a range of thresholds. It has a range 
of 0-1, with a value greater than 0.5 indicating a better-than-random performance event (Fielding, Bell, 1997). 
A rough classification guide is the traditional academic point system (Swets, 1988): poor (0.5-0.6), fair (0.6- 
0.7), good (0.7-0.8), very good (0.8-0.9) and excellent (0.9-1.0). Secondly, the threshold-dependent binomial 
test of omission (Phillips et al., 2006) was employed. The omission test was calculated at a “minimum training 
presence” threshold. At this threshold, the fractional predicted area shows the fraction of all the pixels that were 
predicted suitable for the species. P-values of omission rates in optimal models are considered to be less than 
0.05 (Anderson et al., 2003). 

In addition, for comparative purposes, niche overlap was measured using the Schoener’s index of niche 
breadth (D) (Warren et al., 2008). Indices may range from 0 (indicating that niche models are completely 
different) to 1 (indicating that niche models are identical). 


Results 


Model validation and influencing bioclimatic variables 


The most important sets of variables assessed in performance of models were identified 
separately for the bioclimatic parameters from the WorldClim database and the ENVIREM 
dataset (table 2). 

Maxent produces two ROC curves based on either the initial (training) data and on the 
validation (test) data with associated AUC values. The indicators were obtained using 30 % 
of the training data as test localities for evaluating the performance statistics. The omission 
test was calculated at a Minimum training presence threshold value. At this threshold, the 
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Table 2. Results of Maxent variable selection 





Bioclimatic parameters from the WorldClim database 
BIO 3: Isothermality 

BIO 6: Min Temperature of Coldest Month 

BIO 7: Temperature Annual Range 

BIO 14: Precipitation of Driest Month 

BIO 15: Precipitation Seasonality 

BIO 19: Precipitation of Coldest Quarter 








Species B Variables 
A. farris 1 BIO 3, BIO 7, BIO 15, BIO 19 
A. siro 1 BIO 3, BIO 7, BIO 15, BIO 19 
L. destructor 1 BIO 3, BIO 6, BIO 15, BIO 19 
Gl. domesticus 2 BIO 3, BIO 15, BIO 19 
T. molitor 1 BIO 3, BIO 14, BIO 15 


T. putrescentiae 1 BIO 3, BIO 14 


Climatic variables from the ENVIREM dataset (abbreviations in parentheses) : 
Thornthwaite aridity index (aridityIndex Thornthwaite), 

Climatic moisture index (climaticMoistureIndex), 

Continentality (continentality), 

Emberger’s pluviothermic quotient (embergerQ), 

Minimum temperature of the warmest month (minTempWarmestMonth), 
Maximum temperature of the coldest month (maxTempColdestMonth), 
Monthly variability in potential evapotranspiration (PET)(PETseasonality), 
Mean monthly PET of coldest quarter (PETColdestQuarter) 











Species B Variables 

A. farris 1 aridityIndexThornthwaite, climaticMoistureIndex, continentality, embergerQ, minT- 
empWarmestMonth 

A. siro 1 aridityIndex Thornthwaite, continentality, PETColdestQuarter 

L. destructor 1 aridityIndex Thornthwaite, continentality, PETColdestQuarter, PETseasonality 

Gl. domesticus 1 aridityIndexThornthwaite, continentality, PETColdestQuarter, PETseasonality, max- 
TempColdestMonth 

T. molitor 2 aridityIndexThornthwaite, climaticMoistureIndex, continentality 

T. putrescentiae 1 aridityIndexThornthwaite, embergerQ, PETColdestQuarter, PET seasonality 





fractional predicted area shows the fraction of all the pixels that were predicted suitable for 
the species. Results of threshold-dependent omission test and threshold-independent ROC 
test are presented in table 3 and 3a (statistically significant p-values < 0.05 are marked with 
an asterisk’). 

In terms of the threshold-dependent tests, models incorporating bioclimatic parameters 
from the WorldClim database performed poorly (in all cases p-values above 0.05), whereas 
models using climatic variables from the ENVIREM dataset showed better performance 
(p-values less than 0.05) regarding such species as A. siro, L. destructor and Gl. domesticus. 
The calculated ROC showed that the AUC values of training datasets and test data sets 
were on average higher in models using climatic variables from the ENVIREM dataset, 
particularly for the models considering the mentioned species. These AUC values appear 
to be centred around 0.8, meaning good (0.7-0.8) and/or very good (0.8-0.9) performance 
and indicating the ability of the models to recognize and order a certain point in the 
considered area as suitable or unsuitable. This result was obtained for both the training and 
the test data, with relatively small differences in AUC values (AUC œ suggesting a robust 
performance of the Maxent algorithm to capture the changes in certain environmental 
variables over point localities. 

Table 4 gives estimates of the relative contributions of environmental variables from 
the ENVIREM dataset to the Maxent models for A. siro, L. destructor and Gl. domesticus 
based on permutation importance. 

Based on the jackknifing analysis results, the highest gain in the A. siro model when 
used in isolation is PETcoldQ, which therefore appears to have the most useful information 
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Table 3. Results of threshold-dependent omission test and threshold-independent ROC test 





Models using bioclimatic parameters from the WorldClim database 


























Threshold-dependent test Threshold-independent test 
Species Fractional Test omission | P-value | Test AUC Training AUC | AUC, 
predicted area rate 
A. farris 0.5449 0.0900 0.361 0.6449 0.7249 0.0800 
A. siro 0.8484 0.0820 0.350 0.5850 0.7009 0.1159 
L. destructor 0.8859 0.0320 0.416 0.6369 0.6919 0.0550 
Gl. domesticus 0.8933 0.0267 0.431 0.6499 0.6756 0.0257 
T. molitor 0.6776 0.1000 0.387 0.6982 0.7570 0.0588 
T. putrescentiae 0.7612 0.0300 0.397 0.6526 0.6869 0.0343 





Table 3a. Results of threshold-dependent omission test and threshold-independent ROC test 





Models using climatic variables from the ENVIREM dataset 


























Threshold-dependent test Threshold-independent test 
Species Fractional Test omission | P-value | Test AUC | Training AUC | AUC ae 
predicted area rate 

A. farris 0.5627 0.1133 0.364 0.6206 0.6960 0.0754 
A. siro 0.5852 0.0760 0.048* 0.7545 0.8256 0.0711 
L. destructor 0.5788 0.0413 0.003* 0.7855 0.8304 0.0449 
Gl. domesticus 0.5623 0.0450 0.004* 0.7950 0.8374 0.0424 
T. molitor 0.8700 0.0600 0.423 0.5454 0.5831 0.0377 
T. putrescentiae 0.6504 0.0700 0.330 0.6783 0.7292 0.0509 





by itself. The environmental variable that decreases the gain the most when it is omitted 
is PETcoldQ, which therefore appears to have the most information that is not present in 
the other variables. For the L. destructor model the environmental variable with the highest 
gain when used in isolation is aridityIndexThornthwaite, which therefore appears to have 
the most useful information by itself, whereas the environmental variable that decreases 
the gain the most when it is omitted is PETcoldQ, which therefore appears to have the 
most information that is not present in the other variables. Finally, the environmental 
variable with the highest gain in the Gl. domesticus model when used in isolation is 
aridityIndexThornthwaite, which therefore appears to have the most useful information by 
itself, and the environmental variable that decreases the gain the most when it is omitted 
is PETcoldQ, which therefore appears to have the most information that is not present in 
the other variables. 

In summary, based on the jackknifing analysis and permutation importance, 
the environmental variables that most strongly influenced the habitat suitability for 
the considered mite species are PETcoldQ and the aridityIndexThornthwaite. Their 
joint contribution amounted 88.4 %, 49.4 %, and 42.4 % in the models for A. siro, 
L. destructor and Gl. domesticus, respectively. In addition, continentality was an 
important variable affecting the bioclimatic niche of each of the analyzed species, 
particularly L. destructor. 


Table 4. Contributions of environmental variables (expressed as %) to Maxent model performance 


Variables L. destructor Gl. domesticus 


aridityIndexThornthwaite 41.1 21.1 16.0 
continentality 11.6 24.6 14.0 
maxTempColdestMonth - - 21.7 
PETseasonality - 26.1 21.9 


PETColdestQuarter 47.3 28.3 26.4 
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This strong association is evident in the Maxent response curves. These curves allow 
to see how each environmental variable affects the Maxent prediction by showing how 
the logistic prediction changes as each environmental variable is varied, keeping all other 
environmental variables at their average sample value (figs 4-9). The 10 percentile training 
presence logistic threshold value, which provides a better ecologically significant result 
when compared with more restricted thresholds values (Phillips, Dudik, 2008), was used 
for the cut-off. 


Discussion 


According to the measures of fit applied to determine the ability of niche models 
to predict the distribution of the considered mite species, ENVIREM-based models, at 
least for A. siro, L. destructor and GI. domesticus, outperformed models using bioclimatic 
variables from WorldClim (Hijmans et al., 2005). WorldClim variables are the most 
commonly employed set of environmental data layers used for this purpose, on account of 
their high resolution, global coverage, and availability, however the biological suitability of 
these bioclimatic variables and other such environmental datasets for niche modeling has 
been questioned. In the recently proposed ENVIREM dataset variables have been selected 
on the grounds that they are likely to have direct relevance to ecological or physiological 
processes determining distributions of many species (Title, Bemmels, 2018). They should 
therefore facilitate ecologically-informed variable selection, and may also result in improved 
model performance using statistical variable-thinning approaches. These variables have 
been largely derived from the same underlying dataset as the bioclimatic variables from 
WorldClim (e. g., both “continentality” and “Temperature annual range” are reflections of 
the range of temperatures between warm and cold season), but are intended to complement 
the existing WorldClim dataset (Hijmans et al. 2005) by introducing novel information not 
captured by this database. Through several case studies, P. O. Title and J. B. Bemmels (2018) 
had shown that the ENVIREM variables improved model performance and were valuable 
additions to the set of variables that are currently widely used in species distribution 
modeling. In our modeling exercise the ENVIREM variables too had substantially 
improved model performance by using novel information captured in the PETcoldQ and 
aridityIndexThornthwaite variables. 

A variety of abiotic factors, such as temperature and precipitation, are consistently 
found to be primary determinants of species distributions at broad scales (Wiens, 
2011), but prominent amongst the former is evapotranspiration (ET), often found to 
be one of the best climatic correlates shaping the ecological niche together with species 
distributions and richness (Currie, 1991; Fisher et al., 2011). This was supported in our 
research where PETcoldQ appeared to play a profound role. ET estimates provide an 
indication of ecologically important aspects of climate linked to energy and water supply 
and can be used as measures of ecological energy regimes. The rate of ET depends on the 
intensity of solar radiation, air temperature, humidity and wind speed, and is estimated 
indirectly from meteorological temperature (Clarke, 2017). Potential evapotranspiration 
(PET) is considered a measure of ambient energy and is often correlated with temperature 
variables (Hawkins et al., 2003). Despite this importance, a quantitative synthesis 
analyzing the contribution of over 400 distinct environmental variables to 2040 Maxent 
species’ distribution models PET is poorly represented: summer PET is accounted for in 
34 articles, whereas winter PET only in three (Bradie, Leung, 2017). Strong hump-shaped 
relationships between the probability of presence and PETcoldQ were observed for all 
the three species of mites for which there were reliable SDMs. Probability of presence 
first increased, and then decreased with increasing PETcoldQ, suggesting that moisture 
could be a limiting factor since it is frozen when energy input is at the lowest during 
the coldest months of the year. On the other hand, increasing energy input at this time 
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of the year may cause the release of excessive moisture. This moisture in combination 
with winter coldness is likely to act synergistically and could create critical conditions. 
Therefore there is an apparent optimum of energy input for the cold season under which 
the considered species are able to survive. 

The more complex climatic indices included in the ENVIREM variables (e. g., 
thermicity, aridity, moisture, Emberger’s pluviothermic quotient, monthly variability in 
potential evapotranspiration) are considered to characterize environmental conditions 
that are more directly physiologically relevant to a given species than simple descriptors 
of climate such as temperature or precipitation alone (Title, Bemmels, 2018). One of these, 
the Thornthwaite aridity index, found to be influential in shaping the abiotic niche of the 
considered species, is an indication of moisture deficit and employs evapotranspiration 
and precipitation factors. In the same way, as for the PETcoldQ, strong hump-shaped 
relationships between the probability of presence and aridityIndexThornthwaite were 
observed for the mites species for which there were reliable SDMs, assuming an optimum 
between “too dry” and “too wet” conditions. Another variable from this cohort, the 
monthly variability in potential evapotranspiration (PETseasonality), has been found to be 
important in shaping the abiotic niche in two of the considered mite species: L. destructor 
and GI. domesticus. 

Other influential variables, namely, continentality (calculated as the average 
temperature of warmest month minus the average temperature of coldest month) and 
the maximum temperature of the coldest month (maxTempColdestMonth), are more 
straightforward in their interpretation: increasing continentality reduces the probability 
of presence (i. e., habitat suitability) (fig. 10, exemplified by L. destructor), whereas 
raising maxTempColdestMonth enhances the probability of presence in the model for 
GI. domesticus (fig. 11). 

In terms of the bioclimatic niche the considered mite species display substantial niche 
overlap. Niches of L. destructor and Gl. domesticus are almost identical (Schoener’s D = 
0.929), what could be assumed by the similarity of the response curves and the predictor 
variables they share in common. The niche of A. siro also considerably overlaps with the 
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Fig. 11. Response of Gl. domesticus to maxTempColdestMonth. 
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niches of the two former species, but to a somewhat lesser extent (a 90 % overlap with the 
niche of L. destructor and an 87.5f % overlap with the niche of GI. domesticus). 

Based on our analysis, the bioclimatic variables that influence suitable areas for the 
mite species analyzed seem mostly to be represented by PET-based indices characterizing 
ambient energy, which are often correlated with temperature variables, moisture regimes, 
and strong fluctuations in temperature reflecting the severity of climate and/or extreme 
weather events. In terms of statistical significance, these relationships have been confirmed 
in models created for the species A. siro, L. destructor and Gl. domesticus, for which there 
has been sufficient point data. Matching relationships seem to exist for A. farris, T. molitor 
and T. putrescentiae, however the available data on these species are perhaps too scarce to 
make any robust conclusions. 

In this study we focused on an attempt to highlight important variables shaping the 
current bioclimatic niche of a selected cohort of mite species infesting grain storages, 
realizing that mite persistence is affected by a combination of physical and biological 
factors which interact in a complex manner (Sinha, Wallace, 1973), including interactions 
regarding indoor and outdoor environments. Though microclimates are supposed 
to modulate the responses to the macroclimate, our coarse scale modeling exercise has 
shown the importance and significance of bioclimatic variables in shaping the niche of the 
considered mite species. 

Traditionally, determining such driving factors would require laborious field 
measurements of the key environmental variables in natural populations. However, 
the recent development of remote sensing technologies and rapidly accumulating 
environmental data derived from geographic information systems (GIS) now provide 
information on the patterns of environmental variation at a variety of scales. In 
combination with species occurrence data, these tools permit detailed description of the 
environmental conditions at georeferenced natural population localities for one or more 
species (Nakazato et al., 2010). 

Importantly, the use of SDMs has allowed to identify the environmental and climatic 
features that characterize the species’ niche (their “bioclimate envelope”). On the other hand, 
the SDMs showed that the predictive performance of the models (calculated using AUC) 
never reached values, which could be attributed to the category “excellent” (AUC > 0.9). This 
occurs when ecological parameters are omitted from the modelling framework and lead to 
an inaccurate or insufficient description of the species’ distribution and niche (Hanspach et 
al., 2010). In our case this could be expected as far as only bioclimatic parameters were used 
and it is reasonable to suggest that factors other than climate limit distributions and shape 
the niches of the considered mite species. Including other biologically relevant parameters 
and non-climate variables can contribute important information to SDMs, especially when 
modelers have specific knowledge about how these variables relate to the species (Bucklin 
et al., 2016). 

In addition, broadly characterized as “cosmopolitan”, the considered mite species can 
be hypothesized to be more generalist than ones with smaller home ranges. Because of this 
they are inherently difficult to model as their ecological requirements are less clear. This is 
in line with results from previous studies, which show that species with wider ranges are 
less well predicted (Segurado, Araujo, 2004). 

Stored grain is a man-made ecosystem undergoing persistent interactions between 
several abiotic and biotic factors, and which more or less is affected by the surrounding 
bioclimatic environment. As noticed before, grain storage itself can be heavily affected by 
climate change (Moses et al., 2015), therefore safe food grain storages by themselves are 
considered as a measure to adapt to the changing global climates (UNEP, 2010). Global 
climate patterns have been shown to change notably and are predicted to continue to 
change in the future. Most of the several dozens of predictive models indicate that average 
temperature can increase by 1.7-5.3 °C, as a result of doubling CO, concentration within 
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next 60-100 years. 2.3 °C is a value most often mentioned what means an increase by 0.3 °C 
a decade (IPCC, 2007). 

Predictions based on the emission scenario A1B (Kriticos et al., 2012) and targeting 
the year 2030 show that specific climate changes will affect our study area that are likely to 
have a profound effect on mite development with increases in numbers and growth rates. 
Winter months are expected to be wetter and warmer. Heightening temperatures during 
the coldest weeks will definitely favour the considered species. On the other hand, summers 
are predicted to be hotter and drier, and could lower the average monthly relative humidity 
(RH), which even today (fig. 3) seems to be critical for the warmest months. For practical 
reasons investigations are needed to see how future management targets will be met, based 
on the latest climate scenarios and projected climate changes. 


Conclusions 


The use of SDMs (specifically Maxent, based on the maximum entropy approach) has 
allowed to identify climatic features that characterize the species’ niche (their “bioclimate 
envelope”) of a cohort of mite species known to be associated with stored grains and other 
foods. 

Using the ENVIREM dataset of bioclimatic variables, performance of the most robust 
models was mostly influenced by: 1) indices based on potential evapotranspiration, which 
characterize ambient energy and are mostly correlated with temperature variables, moisture 
regimes, and 2) strong fluctuations in temperature reflecting the severity of climate and/or 
extreme weather events. 

Most likely because only bioclimatic parameters were employed the performance of 
the SDMs (evaluated by the AUC) was < 0.9, reasonably suggesting that factors other than 
climate limit distributions and shape the niches of the considered mite species. Including 
other biologically relevant parameters and non-climate variables could contribute to the 
better performance of the corresponding SDMs and give new insights into the bioecology 
of the species. 

Although the considered mite species occupy man-made ecosystems, they remain 
more or less affected by the surrounding bioclimatic environment and therefore could be 
subjected to contemporary climate change. In this respect investigations are needed to see 
how this will affect future management targets. 
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